!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   FEMilaro -- Finite Element Method toolkit
!!
!! FEMilaro is free software; you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation; either version 3 of the License, or
!! (at your option) any later version.
!!
!! FEMilaro is distributed in the hope that it will be useful, but
!! WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
!! General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with FEMilaro; If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>

!>\brief
!! 
!! Turbulent flux, dynamic Germano model (scalar viscosity)
!!
!! \n
!!
!! The dissipative fluxes are obtained adding a viscous contribution,
!! computed as in \c mod_viscous_flux, and a turbulent contribution.
!! For the turbulent contribution, we use a dynamic model with
!! coefficients averaged on each element.
!!
!! Concerning the momentum flux, we have
!! \f{displaymath}{
!!  \underline{\underline{\mathcal{F}}}^d_{\underline{U}} =
!!  \underbrace{
!!  \underbrace{-\rho\nu_{sgs}\left[ \nabla\underline{u} + \nabla\underline{u}^T -
!!  \frac{2}{d}\nabla\cdot\underline{u}\,\mathcal{I}\right]}_{
!!  \displaystyle\underline{ \underline{\tau}}^d}
!!  +\underbrace{C_I\rho\Delta^2 \left|\nabla\underline{u} +
!!  \nabla\underline{u}^T\right|^2\,\mathcal{I}}_{
!!  \displaystyle\frac{1}{d}\rho e_K\mathcal{I}} }_{\displaystyle 
!!  \displaystyle\underline{ \underline{\tau}}}
!!  +\,\underline{\underline{\mathcal{F}}}^{visc}_{\underline{U}} 
!! \f}
!! where \f$d\f$ is the spatial dimension, \f$\Delta\f$ is a
!! characteristic scale,
!! \f{displaymath}{
!! \nu_{sgs} = C_D\Delta^2 \left|\nabla\underline{u} +
!!  \nabla\underline{u}^T\right|,
!! \f}
!! and the viscous terms are described in \c mod_viscous_flux. Notice
!! that this choice ensures \f${\rm
!! tr}(\underline{\underline{\tau}}^d)=0\f$, so that the trace of the
!! turbulent stress tensor is \f$\rho e_K\f$, which can be regarded as
!! the subgrid kinetic energy. The two dimensionless constants
!! \f$C_D\f$ and \f$C_I\f$ are determined dynamically. To this
!! end, let us introduce the following notation:
!!
!<----------------------------------------------------------------------
module mod_dynamic_iso

!-----------------------------------------------------------------------

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_messages, only: &
   mod_messages_initialized, &
   error, warning, info

 use mod_grid, only: &
   mod_grid_initialized, &
   t_grid

 use mod_base, only: &
   mod_base_initialized, &
   t_base

 use mod_atm_refstate, only: &
   mod_atm_refstate_initialized, &
   t_atm_refstate
 
 use mod_turb_flux_base, only: &
   mod_turb_flux_base_initialized, &
   c_turb_mod, t_turb_diags

 use mod_dgcomp_testcases, only: &
   mod_dgcomp_testcases_initialized, &
   t_phc, phc, coeff_visc

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_dynamic_iso_constructor, &
   mod_dynamic_iso_destructor,  &
   mod_dynamic_iso_initialized, &
   t_viscous_flux

 private

!-----------------------------------------------------------------------

! Module types and parameters

 ! public members

 !> Viscous flux
 !!
 !! The field \c grad stores the gradients of the following quatities:
 !! <ul>
 !!  <li> temperature (second index equal to 1)
 !!  <li> velocity    (second index from 2 to <tt>1+grid\%d</tt>)
 !!  <li> tracers     (second index \f$\geq\f$<tt>2+grid\%d</tt>)
 !! </ul>
 !! Some private fields are also included for the computation of the
 !! diagnostics.
 type, extends(c_turb_mod) :: t_viscous_flux
  real(wp), private, allocatable :: wg(:)
  real(wp), private, allocatable :: base_p(:,:)
 contains 
  procedure, pass(tm) :: init             => vf_init
  procedure, pass(tm) :: grad_diagnostics => vf_grad_diags
  procedure, pass(tm) :: compute_coeffs   => vf_comp_coeffs
  procedure, pass(tm) :: flux             => vf_flux
 end type t_viscous_flux

! Module variables
 ! public members
 logical, protected ::               &
   mod_dynamic_iso_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_viscous_flux'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_dynamic_iso_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_kinds_initialized.eqv..false.) .or. &
    (mod_messages_initialized.eqv..false.) .or. &
        (mod_grid_initialized.eqv..false.) .or. &
        (mod_base_initialized.eqv..false.) .or. &
(mod_atm_refstate_initialized.eqv..false.) .or. &
(mod_turb_flux_base_initialized.eqv..false.).or.&
(mod_dgcomp_testcases_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_dynamic_iso_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   mod_dynamic_iso_initialized = .true.
 end subroutine mod_dynamic_iso_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_dynamic_iso_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_dynamic_iso_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   mod_dynamic_iso_initialized = .false.
 end subroutine mod_dynamic_iso_destructor

!-----------------------------------------------------------------------

 pure subroutine vf_init(tm,td,grid,base,ntrcs)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base
  integer, intent(in) :: ntrcs
  class(t_viscous_flux), intent(inout) :: tm
  type(t_turb_diags), intent(out) :: td

   allocate(tm%wg(base%m)); tm%wg = base%wg
   allocate(tm%base_p(base%pk,base%m)); tm%base_p = base%p

   tm%d           = grid%d
   tm%ntrcs       = ntrcs
   tm%ngrad       = 1 + tm%d + tm%ntrcs
   allocate(tm%grad(tm%d,tm%ngrad,base%pk,grid%ne))
   tm%initialized = .true.

   td%ndiags = 2
   allocate( td%diag_names(td%ndiags) ,  &
     td%diags(td%ndiags,base%pk,grid%ne) )
   td%diag_names(1) = '|nabla^S u + lambda*div(u)|'
   td%diag_names(2) = 'viscous dissipation'
 end subroutine vf_init

!-----------------------------------------------------------------------

 pure subroutine vf_grad_diags(tm,gd,consv,atm_ref)
  class(t_viscous_flux), intent(in) :: tm
  class(t_atm_refstate), intent(in) :: atm_ref(:)
  !> conservative variables (deviations)
  real(wp), intent(in) :: consv(:,:)
  real(wp), intent(out) :: gd(:,:) !< diagnostics

  integer :: id, it
  real(wp), dimension(size(consv,2)) :: rho, e, t
  real(wp), dimension(tm%d,size(consv,2)) :: u, uu

   ! include the reference state
   rho = consv(1,:) + atm_ref%rho
   e   = consv(2,:) + atm_ref%e
   u   = consv(2+1:2+tm%d,:)

   ! velocity
   do id=1,tm%d
     uu(id,:) = u(id,:)/rho
   enddo
   ! temperature
   t = (e/rho - 0.5_wp*sum(uu**2,1) - atm_ref%phi) / phc%cv

   ! collect values
   gd(   1      ,:) = t
   gd(1+1:1+tm%d,:) = uu
   do it=1,tm%ntrcs ! tracers
     gd(1+tm%d+it,:) = consv(1+tm%d+it,:)/rho
   enddo
 
 end subroutine vf_grad_diags

!-----------------------------------------------------------------------

 pure subroutine vf_comp_coeffs(tm,grid,base)
  type(t_grid), intent(in) :: grid
  type(t_base), intent(in) :: base
  class(t_viscous_flux), intent(inout) :: tm
   !> Nothing to do
 end subroutine vf_comp_coeffs

!-----------------------------------------------------------------------

 !> \warning If the optional argument \c td is present, it is assumed
 !! that the function is called for all the quadrature points in one
 !! element (i.e. not for a side).
 pure subroutine vf_flux(fem, tm,ie,x,bp, consv,atm_ref, rho,p,uu,cc, &
                         td)
  class(t_viscous_flux), intent(in) :: tm
  integer,  intent(in) :: ie
  real(wp), intent(in) :: x(:,:)  !< coordinates
  real(wp), intent(in) :: bp(:,:) !< basis functions
  !> conservative variables (deviations)
  real(wp), intent(in) :: consv(:,:)
  class(t_atm_refstate), intent(in) :: atm_ref(:)
  !> Total values for density, pressure, velocty and tracers
  real(wp), intent(in) :: rho(:), p(:), uu(:,:), cc(:,:)
  real(wp), intent(out) :: fem(:,:,:)
  type(t_turb_diags), intent(inout), optional :: td

  integer :: i, l, id, jd
  real(wp) :: &
   nu(2,size(x,2)), & !< viscosity
   gradg(tm%d,tm%ngrad,size(x,2)), & !< gradient in physical space
   gsuu (tm%d,tm%d,    size(x,2)), & !< velocity symmetric gradient
   ldivu, sym(tm%d,tm%d), &
   s_abs(size(x,2)), diss(size(x,2))
 
   ! problem coefficients
   nu = coeff_visc( x , rho,p,p/(rho*phc%rgas) )

   ! gradients
   !do i=1,tm%ngrad
   !  gradg(:,i,:) = matmul(tm%grad(:,i,:,ie),bp)
   !enddo
   gradg = 0.0_wp ! equivalent implementation
   do l=1,size(bp,2)
     do i=1,size(bp,1)
       gradg(:,:,l) = gradg(:,:,l) + tm%grad(:,:,i,ie)*bp(i,l)
     enddo
   enddo
   do id=1,tm%d
     do jd=1,id-1; gsuu(id,jd,:) = gsuu(jd,id,:); enddo ! symmetry
     do jd=id,tm%d
       gsuu(id,jd,:) = gradg(jd,1+id,:)+gradg(id,1+jd,:)
     enddo
   enddo

   do l=1,size(x,2)
     ldivu = (-2.0_wp/3.0_wp)*0.5_wp*tr(gsuu(:,:,l)) ! lambda*div(u)
     sym = gsuu(:,:,l) ! symmetric gradient
     do id=1,tm%d
       sym(id,id) = sym(id,id) + ldivu
     enddo

     ! momentum flux
     fem(:,2:1+tm%d,l) = -rho(l)*nu(2,l)*sym

     ! energy flux
     fem(:,   1    ,l) = -rho(l)*nu(1,l)*phc%cp*gradg(:,1,l)  &
                        + matmul( uu(:,l) , fem(:,2:1+tm%d,l) )

     ! tracers
     fem(:,2+tm%d: ,l) = -rho(l)*nu(1,l)*gradg(:,2+tm%d:,l)
     
     if(present(td)) then
       s_abs(l) = sqrt(sum(sym**2))
       diss (l) = -sum(fem(:,2:1+tm%d,l)*sym) ! use sym^T=sym
     endif
   enddo

   ! the projections are simple because the basis is orthonormal
   if(present(td)) then
     td%diags(1,:,ie) = matmul( tm%base_p , tm%wg*s_abs )
     td%diags(2,:,ie) = matmul( tm%base_p , tm%wg*diss  )
   endif

 end subroutine vf_flux
 
!-----------------------------------------------------------------------

 pure function tr(a)
  real(wp), intent(in) :: a(:,:)
  real(wp) :: tr
  integer :: i

   tr = a(1,1)
   do i=2,minval(shape(a))
     tr = tr + a(i,i)
   enddo
 end function tr

!-----------------------------------------------------------------------

end module mod_dynamic_iso

